function [initial_new,value,iter_2s,diff,nfeval] = Full_3Stages(initial,DP,I_N,A,R,TG,MON,rp,snPass)

iter_2s = 1;
dmc     = 0.05;
diff    = 10;
nfeval  = 0;
MaxIter = 1000;

initialtr = coefftran(initial,1);

if rp.cc~=Inf
    
    while diff>=rp.cc
        
        params=initialtr;
        
        initial1 = initialtr(1:3);
        initial2 = initialtr(4:5);
        initial3 = initialtr(6:8);
        
        %*******************************************************;
        % Stage 1: muT, sdT, theta                              ;
        %*******************************************************;
        stage=1;
        fprintf('\n Iteration %d, Stage %d \n',iter_2s,stage)
        
        options= ...
            optimset('Display','iter','UseParallel','always','MaxIter',MaxIter,'MaxFunEvals',10000,'TolFun',1e-10,'DiffMinChange',0.1);
        
        [theta_hat1,value,flag,info] = fminunc('Full_Likelihood',initial1,options,DP,I_N,A,R,TG,MON,params,rp,snPass,stage);
        
        display('LL & Parameters of Stage 1: muT sigmaT rho')
        [value,theta_hat1]
        nfeval = nfeval + info.funcCount;
        
        params=[theta_hat1,initial2,initial3];
        
        %*******************************************************;
        % Stage 2: muF, fshift, ftreat, fmoni                   ;
        %*******************************************************;
        stage=2;
        fprintf('\n Iteration %d, Stage %d \n',iter_2s,stage)
        
        options= ...
            optimset('Display','iter','UseParallel','always','MaxIter',MaxIter,'MaxFunEvals',10000,'TolFun',1e-10,'DiffMinChange',0.5);
        
        [theta_hat2,value,flag,info] = fminunc('Full_Likelihood',initial2,options,DP,I_N,A,R,TG,MON,params,rp,snPass,stage);
        
        display('LL & Parameters of Stage 2: sF0 sF1')
        [value,theta_hat2]
        nfeval = nfeval + info.funcCount;
        
        params=[theta_hat1,theta_hat2,initial3];
        
        %*******************************************************;
        % Stage 3: sdF0, sdF1                                   ;
        %*******************************************************;
        stage=3;
        fprintf('\n Iteration %d, Stage %d \n',iter_2s,stage)
        
        options= ...
            optimset('Display','iter','UseParallel','always','MaxIter',MaxIter,'MaxFunEvals',10000,'TolFun',1e-10,'DiffMinChange',2);
        
        [theta_hat3,value,flag,info] = fminunc('Full_Likelihood',initial3,options,DP,I_N,A,R,TG,MON,params,rp,snPass,stage);
        
        display('LL & Parameters of Stage 3: meanF params')
        [value,theta_hat3]
        nfeval = nfeval + info.funcCount;
        
        %*******************************************************;
        % Updating initial values and calculating difference for convergence criterion;
        %*******************************************************;
        
        initial_new=[theta_hat1,theta_hat2,theta_hat3];
        
        diff=(initial_new-initialtr)*(initial_new-initialtr)'
        initialtr=initial_new;
        iter_2s=iter_2s+1
        dmc=dmc+0.1;
        
    end
else
    options= ...
        optimset('Display','iter','MaxIter',MaxIter,'MaxFunEvals',10000,'TolFun',1e-10,'DiffMinChange',0.05);
    
    params = [];
    stage  = [];
    
    ub1 = [Inf   50 1  Inf Inf Inf  Inf  Inf];
    lb1 = [-Inf  0  -1 0   0   -Inf -Inf -Inf];
    [initial_new,value,flag,info] = fmincon('Full_Likelihood',initialtr,[],[],[],[],lb1,ub1,[],...
        options,DP,I_N,A,R,TG,MON,params,rp,snPass,stage);
end